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Many  natural  processes  can  be  characterized  by  their  scale-invariance  property. 
In  this  study,  we  present  the  results  of  potential  multiple  scalings  in  the  long- 
term heart  rate  data  from  young  healthy  adults  subjected  to  normal  daily  activity. 
Our  approach  is  based  on  the  direct  check  of  the  probabilistic  structure  of  the 
increment  process.  Results  from  fractional  Brownian  motion  are  compared  and 
the  generating  mechanism  for  multiple  scaling  is  discussed  in  the  context  of  scale- 
invariance  formalism. 
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1 Introduction 

Many  physical  processes  are  scale  invariant,  giving  rise  to  renormalizable  structure 
and  self-similarity  in  space  and  time.  It  is  quite  often  that  the  renormalization 
can  be  conducted  at  different  scales  with  different  scaling  exponents,  namely,  a 
multiple  scaling  property.  Multiple  scaling  was  observed  in  many  natural  phenom- 
ena ranging  from  physical1,2’3,4’5’6 , biological7 ,8’9,10 , to  economical  systems1’11’12. 
When  a multiplicative  mechanism  is  involved  in  the  dynamics,  multifractal  theory 
can  be  used  to  characterize  a (continuous)  set  of  singularity  exponents  and  the 
(Hausdorff)  dimension  of  their  supports  via  a Legendre  transformation.  Successful 
applications  of  the  formalism  revealed  deep  insights  about  many  natural  processes 
such  as  the  density  profile  in  diffusion  limited  aggregates13’14,  the  velocity  and  dis- 
sipation fields  in  fully  developed  turbulence15,2’3,  the  money  exchange  index  from 
financial  market1,11’12  and  network  traffic5’6.  In  general,  there  is  not  an  unified 
theory  for  multiple  scaling  and,  sometimes,  only  finite  number  of  scaling  exponents 
can  be  ascertained.  In  this  work,  we  will  present  such  an  example  in  the  long-term 
heart  rate  variability  (HRV)  from  healthy  young  adults  and  discuss  its  generating 
mechanism. 

The  study  of  HRV  has  continued  to  draw  great  interests  in  recent  years.  It  is 
of  both  fundamental  and  clinical  importance.  In  particular,  a “healthy”  heart  typ- 
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ically  shows  a power-law  like  power  spectrum16 , which  can  imply  scale  invariance0 
and  “self-similarity”  in  the  autonomous  nervous  system.  Losing  such  variability 
was  found  to  correlate  well  with  the  mortality  rate  of  heart  diseased  patients17. 
The  variation  in  the  HRV  scaling  has  been  discussed  by  a number  of  researchers 
in  the  past7,8,10.  In  particular,  Peng  et.  al.  developed  the  systematic  method  and 
compared  the  power-law  scalings  in  the  very  short  time  scale  (a  few  heart  beats) 
with  the  “asymptotic”  behaviour  (above  1000  heart  beats).  Significant  difference 
was  concluded  in  these  widely  separated  time  scales  and  the  characteristic  was 
found  sufficient  to  distinguish  between  the  healthy  subjects  and  patients  with  con- 
gestive heart  failure10.  Di  Rienzo  et.  al.  also  noted  the  variation  of  the  power-law 
exponent  of  HRV  in  frequency,  which  in  turns  causes  variations  of  the  power-law 
scaling  of  the  blood  pressure  as  well7.  Later,  Viswanathan  et.  al.  tackled  the  non- 
uniformity of  power  law  scaling  between  the  healthy  subjects  and  those  with  severe 
heart  disease  from  a point-process  aspect18.  Again,  the  conclusion  was  drawn  on 
the  distinguishability  between  the  health  and  disease.  In  these  previous  studies, 
multiple  scaling  within  physiologically  relevant  range  of  healthy  humans  has  not 
been  explored. 

In  this  work,  multiple  scaling  was  studied  based  on  the  data  increment  process. 
The  increment  process  at  various  time  lags  enables  us  to  focus  on  the  local  scale 
invariance  property  of  the  data.  The  numerical  method  was  developed  and  tested 
on  artificial  time  series  with  one-  and  multiple-  scaling  exponents  and  then  applied 
to  the  heart  rate  data  set.  We  also  compared  with  the  existing  methodology  in 
the  literature  and  found  noticeable  difference  in  the  result.  This  paper  is  organized 
as  follows.  The  main  ideas  and  the  numerical  method  are  introduced  in  the  next 
section.  The  results  of  artificial  time  series  and  the  evidence  of  multiple  scaling 
in  heart  rate  data  are  presented  in  section  3.  In  the  last  section,  the  mechanism 
generating  the  multiple  scaling  in  HRV  will  be  discussed  in  the  context  of  scale 
invariance  formalism. 

2 Extract  Scale  Invariance  from  the  Probability  Density  Function 

In  this  section,  we  will  first  recall  the  basic  definition  of  scale  invariance  of  a pro- 
cess and  then  describe  the  numerical  method  to  extract  multiple  scaling.  Given  a 
time  series,  r(s),  its  scale  invariance  is  defined  by  the  family  of  probability  density 
functions  of  the  increment,  Ar (t;s)  = r(t  + s)  — r(s):  for  any  A and  a fixed  constant 
h,  one  has 

fxt(n)  = \-hft(\~hn)  (1) 

where  ft  denotes  the  probability  density  function  of  n = Ar (t;  s)  and  the  constant 
h is  the  scaling  index  or  scaling  exponent  of  r(s). 

To  study  (1)  numerically  requires  the  assumption  of  stationarity.  This  appears 
to  hold  in  our  application.  When  the  increment  is  stationary,  /t(n)  can  then  be 
estimated  from  the  ensemble  of  {Ar(f;s),s  = 1 , 2,  - ■ . Due  to  the  data  fluctua- 

tion, it  was  found  more  efficient  to  fit  ft(n)  with  a specific  form,  say  gt(n),  than 

“Note  that  power-law  spectrum  is  only  a necessary  condition19,20. 
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Figure  1.  (a)  The  set  Af(T;  h)  plotted  on  the  h — n plane  for  T = {2k , k = 0,1, 2, 3}  and 
h0  = 0.5,  dark  pixels:  |A7|  = 0,  white  pixels:  |JVj  > 0,  (b)  fj,(T\h)  vs.  h plot.  The  family 
of  density  functions  used  in  this  demonstration  was  ideal  in  that  they  have  the  Gaussian  form: 
{3t(n)  = exp(— n2/2<r(t)2)/ \J2k cr(i)2,  <r(i)  = \ft\  t = 2k , k = 0,  • • • 11}. 


working  directly  with  the  histogram  of  Ar (t)  in  (1).  For  example,  for  the  fractional 
Brownian  motion,  a Gaussian  form  is  assumed  for  gt(n).  In  this  case,  the  maximum 
likelihood  method  was  used  to  minimize  the  likelihood  function  — log(II;gi(ni))  in 
order  to  extract  the  mean  and  variance  of  a Gaussian  probability  density  function. 
This  approach  is  subjected  to  less  bias  comparing  to,  say,  minimizing  the  L 2 norm 
between  the  histogram  and  gt(n)21.  After  gt{n)  is  defined,  (1)  can  be  studied  by 
systematically  varying  n,  A and  h0  values.  Let 


u(n,  A;  h) 


f\t(n) 


(2) 


and  denote  the  estimated  slope,  dlog(u,)/diog(A),  as  h! . r is  said  to  be  renormaliz- 
able  or  is  scale  invariant  with  a scaling  exponent  h„  if  h!  ~ h0. 

The  parameters  used  in  the  numerical  experiment  were  t = 1,  A £ {2 k,k  = 
0,  • • ■ , 11}  and  h £ {0.01  • k,  k — 0,  • • • , 100}.  The  estimated  slope  b!  is  in  general 
a function  of  n,  h,  and  A.  The  range  of  A in  which  h!  is  estimated  reveals  the 
time  scale  of  the  local  scale  invariance  property  of  r(s).  Once  the  parameters  are 
defined,  we  first  construct,  for  a given  h,  AT(T;  h)  = {n,  \h! {n , A;  h)  — h | < e}  with 
e = 0.01  to  keep  the  error  between  h'  and  h less  less  than  1%.  Here,  T denotes 
the  set  of  A’s  where  the  condition  \h!  — h\  < e is  satisfied.  For  any  given  h and 
T,  A f contains  the  set  of  n’s  in  which  local  scale  invariance  is  defined.  The  scaling 
interval  of  h!  can  thus  be  calculated  by  the  logarithm  of  the  ratio  of  the  largest  and 
smallest  A in  T,  For  each  h,  (1)  also  contains  an  isolated  solution  (when  the  set 
A/"  has  only  one  element)  which  has  no  relevance  in  the  current  context  (Fig.  la). 
Fortunately,  these  isolated  solutions  form  a small  set.  Hence,  given  any  T and  h, 
/.t(T;  h)  = |A/"(T;  h)\,  where  | • | returns  the  number  of  elements  at  a specific  h value, 
will  show  a peak  at  the  desired  scaling  exponent  (Fig.  lb). 
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Figure  2.  (a)  Typical  RR-interval  record  (r (s)  (sec)  vs.  a (xlO3  beat)),  (b)  the  power  spectral 
density  function  of  r(s)  (in  log-log  scale),  (c)  histogram  of  Ar(t),t  = 4 (in  linear-log  scale). 


3 Numerical  Evidence  of  Multiple  Scaling  in  HRV 

3.1  Experimental  Procedure 

Six  young  adults  (average  age:  25  yr,  height:  173  cm,  weight:  74  kg)  participated  in 
the  experiment.  The  subjects  were  allowed  to  conduct  their  normal  daily  activity. 
The  difference  in  the  body  surface  potential  was  sampled  at  1000  Hz  for  a period  of 
approximately  24  hours.  The  data  was  then  down-loaded  to  a PC  and  a specialized 
computer  code  was  written  to  search  for  the  QRS  complex  for  each  heart  beat 
with  proper  filtering  for  events  such  as  skip-beat,  PVC,  and  so  on.  The  time  span 
between  the  successive  contractions  of  the  ventricles,  measured  as  the  RR-interval 
(RRi),  was  finally  extracted  and  used  in  the  scale  invariance  analysis. 

A representative  day-time  RRi  data  is  shown  in  Fig.  2a.  Scale  invariance  in  the 
RRi  data  may  be  implied  by  the  power-law  trend  of  the  power  spectrum  (Fig.  2b). 
But  it  is  clear  that  a single  power-law  is  not  sufficient  to  describe  the  spectrum  over 
the  full  frequency  range. 

3.2  Numerical  Results 

Before  applying  the  numerical  method  to  the  RRi  data,  we  first  tested  it  on  the 
fractional  Brownian  motion  of  scaling  exponent  h0  = 0.26,  0.5,  0.78.  Each  artificial 
time  series  had  40,000  samples  and  was  generated  by  using  the  spectral  method19,20 . 
Fig.  3a  shows  a typical  case  of  the  set  Af  plotted  on  the  h — T plane  where  T = 
2k,k  = 1,  2,  • • ■ Sometimes,  it  is  useful  to  consider  A*(T;  h)  over  all  T sets  for  a 
given  h.  The  result  is  shown  in  Fig.  3b  where  the  scaling  exponent  is  indicated  at 
the  location  of  the  peak.  It  is  clear  that  the  desired  scaling  exponents  were  captured 
with  good  accuracy. 

We  next  apply  the  spectral  method  to  construct  a two-exponent  artificial  time 
series  of  216  data  points  (Fig.  4a).  The  local  scale  invariance  in  the  time  series 
was  defined  by  the  scaling  exponent  h0  = 0.26  for  frequency  below  0. 0025Hz  and 
ha  = 0.78  above  0. 0025Hz.  Fig.  4d  was  plotted  with  ^p(T;/i)  versus  h and  the 
estimated  h's  were  within  3%  of  the  exact  values.  In  Fig.  4c,  we  checked  the 
minimum  time  scale  of  the  exponents,  which  is  defined  by  the  smallest  element  of 
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Figure  3.  (a)  The  set  y,  on  the  h — T plane,  dark  pixels:  |/i|  = 0,  white  pixels:  |/^|  > 0 for  h0  — 0.78. 
(b)  ^(T;  h)  vs.  h for  h0  = 0.26,  0.5,  0.78. 


the  T set.  This  is  obtained  by  constructing  the  set  rj{T\h)  where  T is  the  smallest 
A in  the  corresponding  T set.  The  set  r)  shown  on  the  h — log2(  T)  plane  reveals  the 
minimum  time  scale  from  which  the  exponent  is  defined  (Fig.  4c).  For  the  artificial 
time  series,  it  was  found  that  correct  time  scales  of  the  exponent  were  captured: 
i.e.,  the  “fast”  dynamics  (above  0. 0025Hz)  of  hQ  — 0.78  prevails  in  small  A’s  and 
the  “slow”  dynamics  of  h0  = 0.26  at  large  A (see  Fig.  4c  caption). 

To  study  the  RIti  data,  the  following  form  of  probability  density  function  was 
assumed  (based  on  Fig.  lc): 

_ , , exp(— |n|7^) 

flt(n)  = . 

7 

The  RRi  data  of  all  the  test  subjects  exhibit  multiple  scaling  characteristics.  In 
what  follows,  we  will  present  the  evidence  from  a typical  individual  whose  data  has 
been  given  in  Fig.  2.  In  Fig.  5a,  multiple  peaks  were  seen  from  the  vs.  h 
plot,  indicating  multiple  scaling  of  HRV.  At  least  five  “significant”  exponents  were 
identified  (a  ~ e in  Fig.  5a).  The  minimum  time  scale  of  these  exponents  scattered 
over  the  range  of  A = 21  to  28  (Fig.  5b).  The  inverse  of  this  range  covers  the 
“frequency”  < 0.0039  (1/beat)  to  0.5  (1/beat)  in  the  power  spectrum  (Fig.  2b).  In 
this  range,  two  linear  regions  of  slopes  £ ~ —1.4  and  -2.1  may  be  roughly  defined. 
Based  on  £ — 1 + 2 h,  they  match  nicely  to  two  of  the  peak  locations  h = 0.17  and 
0.6  in  Fig.  5a.  The  size  of  the  scaling  interval  is  shown  by  the  set  Af  plotted  on  the 
h — | T | plane  (Fig.  5c).  A wide  range  of  scaling  interval  associated  with  the  scaling 
exponents  found  in  Fig.  5a  is  again  seen. 

Finally,  we  compared  our  method  with  the  detrend  fluctuation  analysis 
(DFA)8’10  and  found  noticeable  difference.  In  Fig.  6 shows  the  double-logarithmic 
plot  of  the  average  variance  versus  the  scaled  window  length.  Only  one  exponent 
can  be  ascertained  here  since  the  local  scale  invariance  characteristics  have  been 
averaged  out  in  the  process  of  DFA. 


Figure  4.  (a)  Two-time-scale  artificial  time  series,  (b)  power  spectral  density  of  the  time  series: 
the  scaling  exponents  are  h0  = 0.26  for  / < 0.0025ffz  and  h0  = 0.78  for  / > 0.0025-ffz.  (Note: 
0.0025  = 2~8'6).  The  line  segments  have  slops  of  -1.52  and  -2.56,  respectively,  (c)  The  set  r) 
plotted  on  h — log2(T)  plane,  dark  pixels:  |j;j  = 0,  white  pixels:  |r;|  > 0.  The  scaling  range  for 
h = 0.77  was  captured  for  T < 27  and  that  for  h = 0.27  for  T — 27  ~ 29.  (d)  p( T;h)  vs.  h 


Figure  5.  (a)  ^/r(T -,h)  vs.  h plot.  Identified  peaks  at  h=0.17,  0.24,  0.31,  0.44,  0.6  are  labeled  as 
a,  b,  c,  d.  e , respectively;  (b)  the  set  r]  plotted  on  h — log(  T)  plane,  dark  pixels:  |»j|  = 0 and  white 
pixels:  [77]  > 0,  (c)  The  set  of  A7  plotted  on  h — log2(T)  plane,  dark  pixels:  |A/"|  = 0,  white  pixels: 

w\  > 0. 


4 Discussion  and  Future  Outlook 


Numerical  evidence  of  multiple  scaling  of  HRV  in  physiologically  relevant  range 
has  been  presented  for  the  case  of  healthy  young  adults.  The  scale  invariance  was 
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Figure  6.  Detrend  fluctuation  analysis  on  the  RRi  data  shown  in  Fig.  2a:  averaged  variance  of 
the  detrend  time  series  vs.  window  length  plot  on  double  logarithmic  scales  (symbol:  ’o’).  The 
fitted  line  of  slope  ~ 1.05,  suggesting  h = 0.05,  was  shown  as  the  solid  line. 


extracted  based  on  the  property  of  the  family  of  RRi  increment  probability  density 
functions.  This  approach  enables  us  to  explore  local  scale  invariance  property  in 
the  data.  Although  it  is  only  a bit  costlier  in  computation,  details  previously 
unavailable  from  other  global  methods,  such  as  the  power  spectral  density  function 
or  DFA,  may  be  obtained.  It  should  be  noted  that  the  current  methodology  is 
not  able  to  extract  temporal  structure  of  very  short  duration.  When  such  a rare 
event  occurs,  its  characteristics  will  yield  to  that  generated  by  the  more  “regular” 
dynamics  in  the  process  of  estimating  the  density  function.  What  this  work  has 
shown  is  that  even  the  “regular”  neuro-control  of  our  cardiovascular  system  is  rich 
enough  to  encompass  a wide  range  of  time  scales  and  scaling  structures,  i.e. , multiple 
scaling.  Indeed,  the  process  is  very  complex,  as  Hausdorff  and  Peng  showed  from 
artificial  time  series  that  the  “balance”  of  different  inputs  to  the  system  is  also 
crucial  to  generate  the  1/f  scaling8. 

More  detailed  characterization  of  the  density  function  (4)  is  underway  and  will 
provide  the  insight  of  the  generating  mechanism  of  multiple  scaling.  For  example, 
the  variation  of  the  parameters  a,  f3  and  7 in  (4)  can  imply  self-similarity  of  the 
individual  density  function  at  specific  time  lag  (At).  This  in  turns  can  lead  to  scale 
invariant  solution  of  (1)  over  finite  n-interval.  In  particular,  given  a A £ R and 
n'  £ 7n(Af),  assume  the  graph  f\t{n)  is  self-similar  in  the  interval 

AA/xt(AAn')  = fxt{n').  (5) 

Since  it  can  be  written 

f\t(n)  = fxt(Xhn"),  (6) 

(5)  implies 


AA/At(AAn')  = A ~hft(n") 


(7) 
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where  n!  — A hn" . Re-arranging  terms  in  (7),  one  finds 

\h+Afxt(\h+An")  = ft(n")  (8) 

which  is  of  the  same  form  as  (1).  Hence,  (8)  implies  the  existence  of  a new  exponent 
h + A.  Substituting  the  newly  found  exponent  into  the  scale  invariance  formalism 
(1),  even  more  can  be  revealed  by  following  the  same  arguments.  In  general,  we 
found,  when  certain  conditions  are  met,  the  number  of  exponents  can  “grow”  as  a 
power-law  by  repeating  (5)  with  all  the  exponents.  However,  the  scaling  interval 
In  also  shrinks  at  the  same  time  by  a power-law,  making  most  of  the  exponents 
“unobservable22”  . We  will  report  more  details  in  a future  publication. 
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